{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's start with our crazy stock list of imports and setup our environment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# See requirements.txt to set up your dev environment.\n",
    "import os\n",
    "import sys\n",
    "import utm\n",
    "import json\n",
    "import scipy\n",
    "import overpy\n",
    "import urllib\n",
    "import datetime \n",
    "import urllib3\n",
    "import rasterio\n",
    "import subprocess\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from osgeo import gdal\n",
    "from planet import api\n",
    "from planet.api import filters\n",
    "from traitlets import link\n",
    "import rasterio.tools.mask as rio_mask\n",
    "from shapely.geometry import mapping, shape\n",
    "from IPython.display import display, Image, HTML\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.image as mpimg\n",
    "#from scipy import ndimage\n",
    "import warnings\n",
    "from osgeo import gdal\n",
    "\n",
    "\n",
    "from osmapi import OsmApi\n",
    "from geopy.geocoders import Nominatim\n",
    "\n",
    "urllib3.disable_warnings()\n",
    "from ipyleaflet import (\n",
    "    Map,\n",
    "    Marker,\n",
    "    TileLayer, ImageOverlay,\n",
    "    Polyline, Polygon, Rectangle, Circle, CircleMarker,\n",
    "    GeoJSON,\n",
    "    DrawControl\n",
    ")\n",
    "\n",
    "%matplotlib inline\n",
    "# will pick up api_key via environment variable PL_API_KEY\n",
    "# but can be specified using `api_key` named argument\n",
    "api_keys = json.load(open(\"apikeys.json\",'r'))\n",
    "client = api.ClientV1(api_key=api_keys[\"PLANET_API_KEY\"])\n",
    "gdal.UseExceptions()\n",
    "api = overpy.Overpass()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's bring up our slippy map once again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Basemap Mosaic (v1 API)\n",
    "mosaicsSeries = 'global_quarterly_2017q1_mosaic'\n",
    "# Planet tile server base URL (Planet Explorer Mosaics Tiles)\n",
    "mosaicsTilesURL_base = 'https://tiles0.planet.com/experimental/mosaics/planet-tiles/' + mosaicsSeries + '/gmap/{z}/{x}/{y}.png'\n",
    "# Planet tile server url\n",
    "mosaicsTilesURL = mosaicsTilesURL_base + '?api_key=' + api_keys[\"PLANET_API_KEY\"]\n",
    "# Map Settings \n",
    "# Define colors\n",
    "colors = {'blue': \"#009da5\"}\n",
    "# Define initial map center lat/long\n",
    "center = [45.5231, -122.6765]\n",
    "# Define initial map zoom level\n",
    "zoom = 13\n",
    "# Set Map Tiles URL\n",
    "planetMapTiles = TileLayer(url= mosaicsTilesURL)\n",
    "# Create the map\n",
    "m = Map(\n",
    "    center=center, \n",
    "    zoom=zoom,\n",
    "    default_tiles = planetMapTiles # Uncomment to use Planet.com basemap\n",
    ")\n",
    "# Define the draw tool type options\n",
    "polygon = {'shapeOptions': {'color': colors['blue']}}\n",
    "rectangle = {'shapeOptions': {'color': colors['blue']}} \n",
    "\n",
    "# Create the draw controls\n",
    "# @see https://github.com/ellisonbg/ipyleaflet/blob/master/ipyleaflet/leaflet.py#L293\n",
    "dc = DrawControl(\n",
    "    polygon = polygon,\n",
    "    rectangle = rectangle\n",
    ")\n",
    "# Initialize an action counter variable\n",
    "actionCount = 0\n",
    "AOIs = {}\n",
    "\n",
    "# Register the draw controls handler\n",
    "def handle_draw(self, action, geo_json):\n",
    "    # Increment the action counter\n",
    "    global actionCount\n",
    "    actionCount += 1\n",
    "    # Remove the `style` property from the GeoJSON\n",
    "    geo_json['properties'] = {}\n",
    "    # Convert geo_json output to a string and prettify (indent & replace ' with \")\n",
    "    geojsonStr = json.dumps(geo_json, indent=2).replace(\"'\", '\"')\n",
    "    AOIs[actionCount] = json.loads(geojsonStr)\n",
    "    \n",
    "# Attach the draw handler to the draw controls `on_draw` event\n",
    "dc.on_draw(handle_draw)\n",
    "m.add_control(dc)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's review from last time. \n",
    "* We'll query the Planet API and get a list of scenes.\n",
    "* We'll then use pandas and shapely to clean up and filter the results\n",
    "* We'll then render the footprints of the good scenes over our AOI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print AOIs[1]\n",
    "myAOI = AOIs[1][\"geometry\"]\n",
    "\n",
    "# build a query using the AOI and\n",
    "# a cloud_cover filter that excludes 'cloud free' scenes\n",
    "\n",
    "old = datetime.datetime(year=2017,month=1,day=1)\n",
    "\n",
    "query = filters.and_filter(\n",
    "    filters.geom_filter(myAOI),\n",
    "    filters.range_filter('cloud_cover', lt=5),\n",
    "    filters.date_range('acquired', gt=old)\n",
    ")\n",
    "\n",
    "# build a request for only PlanetScope imagery\n",
    "request = filters.build_search_request(\n",
    "    query, item_types=['PSScene3Band']\n",
    ")\n",
    "\n",
    "# if you don't have an API key configured, this will raise an exception\n",
    "result = client.quick_search(request)\n",
    "scenes = []\n",
    "planet_map = {}\n",
    "for item in result.items_iter(limit=500):\n",
    "    planet_map[item['id']]=item\n",
    "    props = item['properties']\n",
    "    props[\"id\"] = item['id']\n",
    "    props[\"geometry\"] = item[\"geometry\"]\n",
    "    props[\"thumbnail\"] = item[\"_links\"][\"thumbnail\"]\n",
    "    scenes.append(props)\n",
    "scenes = pd.DataFrame(data=scenes)\n",
    "# now let's clean up the datetime stuff\n",
    "# make a shapely shape from our aoi\n",
    "portland = shape(myAOI)\n",
    "footprints = []\n",
    "overlaps = []\n",
    "# go through the geometry from our api call, convert to a shape and calculate overlap area.\n",
    "# also save the shape for safe keeping\n",
    "for footprint in scenes[\"geometry\"].tolist():\n",
    "    s = shape(footprint)\n",
    "    footprints.append(s)\n",
    "    overlap = 100.0*(portland.intersection(s).area / portland.area)\n",
    "    overlaps.append(overlap)\n",
    "# take our lists and add them back to our dataframe\n",
    "scenes['overlap'] = pd.Series(overlaps, index=scenes.index)\n",
    "scenes['footprint'] = pd.Series(footprints, index=scenes.index)\n",
    "# now make sure pandas knows about our date/time columns.\n",
    "scenes[\"acquired\"] = pd.to_datetime(scenes[\"acquired\"])\n",
    "scenes[\"published\"] = pd.to_datetime(scenes[\"published\"])\n",
    "scenes[\"updated\"] = pd.to_datetime(scenes[\"updated\"])\n",
    "\n",
    "scenes = scenes[scenes['overlap']>0.9]\n",
    "\n",
    "\n",
    "print len(scenes)\n",
    "# now let's clean up the datetime stuff\n",
    "# make a shapely shape from our aoi\n",
    "portland = shape(myAOI)\n",
    "footprints = []\n",
    "overlaps = []\n",
    "# go through the geometry from our api call, convert to a shape and calculate overlap area.\n",
    "# also save the shape for safe keeping\n",
    "for footprint in scenes[\"geometry\"].tolist():\n",
    "    s = shape(footprint)\n",
    "    footprints.append(s)\n",
    "    overlap = 100.0*(portland.intersection(s).area / portland.area)\n",
    "    overlaps.append(overlap)\n",
    "# take our lists and add them back to our dataframe\n",
    "scenes['overlap'] = pd.Series(overlaps, index=scenes.index)\n",
    "scenes['footprint'] = pd.Series(footprints, index=scenes.index)\n",
    "# now make sure pandas knows about our date/time columns.\n",
    "scenes[\"acquired\"] = pd.to_datetime(scenes[\"acquired\"])\n",
    "scenes[\"published\"] = pd.to_datetime(scenes[\"published\"])\n",
    "scenes[\"updated\"] = pd.to_datetime(scenes[\"updated\"])\n",
    "\n",
    "# first create a list of colors\n",
    "colors = [\"#ff0000\",\"#00ff00\",\"#0000ff\",\"#ffff00\",\"#ff00ff\",\"#00ffff\"]\n",
    "# grab our scenes from the geometry/footprint geojson\n",
    "footprints = scenes[\"geometry\"].tolist()\n",
    "# for each footprint/color combo\n",
    "\n",
    "for footprint,color in zip(footprints,colors):\n",
    "    # create the leaflet object\n",
    "    feat = {'geometry':footprint,\"properties\":{\n",
    "            'style':{'color': color,'fillColor': color,'fillOpacity': 0.1,'weight': 1}},\n",
    "            'type':u\"Feature\"}\n",
    "    # convert to geojson\n",
    "    gjson = GeoJSON(data=feat)\n",
    "    # add it our map\n",
    "    m.add_layer(gjson)\n",
    "# now we will draw our original AOI on top \n",
    "feat = {'geometry':myAOI,\"properties\":{\n",
    "            'style':{'color': \"#FFFFFF\",'fillColor': \"#FFFFFF\",'fillOpacity': 0.1,'weight': 2}},\n",
    "            'type':u\"Feature\"}\n",
    "gjson = GeoJSON(data=feat)\n",
    "m.add_layer(gjson)   \n",
    "m "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now we'll add in our boiler plate activation code for reference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_products(client, scene_id, asset_type='PSScene3Band'):    \n",
    "    \"\"\"\n",
    "    Ask the client to return the available products for a \n",
    "    given scene and asset type. Returns a list of product \n",
    "    strings\n",
    "    \"\"\"\n",
    "    out = client.get_assets_by_id(asset_type,scene_id)\n",
    "    temp = out.get()\n",
    "    return temp.keys()\n",
    "\n",
    "def activate_product(client, scene_id, asset_type=\"PSScene3Band\",product=\"analytic\"):\n",
    "    \"\"\"\n",
    "    Activate a product given a scene, an asset type, and a product.\n",
    "    \n",
    "    On success return the return value of the API call and an activation object\n",
    "    \"\"\"\n",
    "    temp = client.get_assets_by_id(asset_type,scene_id)  \n",
    "    products = temp.get()\n",
    "    if( product in products.keys() ):\n",
    "        return client.activate(products[product]),products[product]\n",
    "    else:\n",
    "        return None \n",
    "\n",
    "def download_and_save(client,product):\n",
    "    \"\"\"\n",
    "    Given a client and a product activation object download the asset. \n",
    "    This will save the tiff file in the local directory and return its \n",
    "    file name. \n",
    "    \"\"\"\n",
    "    out = client.download(product)\n",
    "    fp = out.get_body()\n",
    "    fp.write()\n",
    "    return fp.name\n",
    "\n",
    "def scenes_are_active(scene_list):\n",
    "    \"\"\"\n",
    "    Check if all of the resources in a given list of\n",
    "    scene activation objects is read for downloading.\n",
    "    \"\"\"\n",
    "    retVal = True\n",
    "    for scene in scene_list:\n",
    "        if scene[\"status\"] != \"active\":\n",
    "            print \"{} is not ready.\".format(scene)\n",
    "            return False\n",
    "    return True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now we'll activate our scenes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "to_get = scenes[\"id\"][0:7].tolist()\n",
    "activated = []\n",
    "# for each scene to get\n",
    "for scene in to_get:\n",
    "    # get the product \n",
    "    product_types = get_products(client,scene)\n",
    "    for p in product_types:\n",
    "        # if there is a visual product\n",
    "        if p == \"visual\": # p == \"basic_analytic_dn\"\n",
    "            print \"Activating {0} for scene {1}\".format(p,scene)\n",
    "            # activate the product\n",
    "            _,product = activate_product(client,scene,product=p)\n",
    "            activated.append(product)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# And then download them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tiff_files = []\n",
    "asset_type = \"_3B_Visual\"\n",
    "# check if our scenes have been activated\n",
    "if True:#scenes_are_active(activated):\n",
    "    for to_download,name in zip(activated,to_get):\n",
    "        # create the product name\n",
    "        name = name + asset_type + \".tif\"\n",
    "        # if the product exists locally\n",
    "        if( os.path.isfile(name) ):\n",
    "            # do nothing \n",
    "            print \"We have scene {0} already, skipping...\".format(name)\n",
    "            tiff_files.append(name)\n",
    "        elif to_download[\"status\"] == \"active\":\n",
    "            # otherwise download the product\n",
    "            print \"Downloading {0}....\".format(name)\n",
    "            fname = download_and_save(client,to_download)\n",
    "            tiff_files.append(fname)\n",
    "            print \"Download done.\"\n",
    "        else:\n",
    "            print \"Could not download, still activating\"\n",
    "else:\n",
    "    print \"Scenes aren't ready yet\"\n",
    "\n",
    "sorted(tiff_files)\n",
    "print tiff_files "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's get going with Open Street Maps.\n",
    "* [Open Street Maps](https://www.openstreetmap.org/) is a huge and open collection of data about the Earth. \n",
    "* OSM is free to query. The interfaces are powerful, but hella cryptic. \n",
    "* Let's say we had a pixel in an image and we wanted to know what in the world was at that pixel. \n",
    "* We can use the Open Street Maps [Nominatim](http://wiki.openstreetmap.org/wiki/Nominatim) function to look up what is there, like Google maps.\n",
    "* We can also use the OSM interface to find the 'nodes' near our pixel. \n",
    "* OSM Nominatim works through Lat Long values. To get these lat long values we are going to through [UTM coordinates](https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system).\n",
    "* To get correct the UTM values we'll need to ask GDAL what our UDM zone is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "infile = tiff_files[0]\n",
    "# Open the file\n",
    "gtif = gdal.Open(infile)\n",
    "# Get the project reference object this knows the UTM zone\n",
    "reff = gtif.GetProjectionRef()\n",
    "# arr is the actual image data.\n",
    "arr = gtif.ReadAsArray()\n",
    "# Trans is our geo transfrom array. \n",
    "trans = gtif.GetGeoTransform()\n",
    "# print the ref object\n",
    "print reff\n",
    "# find our UTM zone\n",
    "i = reff.find(\"UTM\")\n",
    "print reff[i:i+12]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now we are going to write a function to convert pixels to UTM\n",
    "* Also a quick function to plot a point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def pixel2utm(ds, x, y):\n",
    "    \"\"\"\n",
    "    Returns utm coordinates from pixel x, y coords\n",
    "    \"\"\"\n",
    "    xoff, a, b, yoff, d, e = ds.GetGeoTransform()\n",
    "    xp = a * x + b * y + xoff\n",
    "    yp = d * x + e * y + yoff\n",
    "    return(xp, yp)\n",
    "\n",
    "def draw_point(x,y,img):\n",
    "    t = 20\n",
    "    # a cloud_cover filter that ex\n",
    "    img[y-t:y+t,x-t:x+t,:] = [255,0,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now let's query a point on our scene and see what OSM tells us.\n",
    "* First we'll define a pixel position\n",
    "* We'll use GDAL to open the scene and then map a pixel to UTM\n",
    "* We'll then convert the UTM value to Lat / Lon using the UTM region we found before. \n",
    "* Then we'll instantiate a Nominatim object and perform a revers lookup and print the results. \n",
    "* We'll then use the OSM Api to get node at this place. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pos = [3000,1400] # this is the pixel we want info abou\n",
    "ds = gdal.Open(infile)\n",
    "# take the GDAL info and make it into UTM\n",
    "my_utm = pixel2utm(ds,pos[0],pos[1])\n",
    "# convert UTM into Lat Long\n",
    "# need to figure out how to get zone info\n",
    "my_lla = utm.to_latlon(my_utm[0],my_utm[1],10,\"N\")\n",
    "# do the lat long look up from OSM\n",
    "geolocator = Nominatim()\n",
    "# reverse look up the are based on lat lon\n",
    "location = geolocator.reverse(\"{0},{1}\".format(my_lla[0],my_lla[1]))\n",
    "# print location info\n",
    "print location.address\n",
    "print location.raw\n",
    "# get the OSM ID info\n",
    "osm_id = int(location.raw[\"place_id\"])\n",
    "print osm_id\n",
    "# create an interface to the OSM API\n",
    "MyApi = OsmApi()\n",
    "# Look up our position \n",
    "print MyApi.NodeGet(osm_id)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now for completeness we'll plot our scene and add the annotation about the spot we found. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.patches import Circle\n",
    "fig,ax = plt.subplots(1)\n",
    "\n",
    "# create our plot\n",
    "plt.imshow(arr[:3,:,:].transpose((1, 2, 0)))#, extent=extent)\n",
    "fig = plt.gcf()\n",
    "# add our annotation\n",
    "plt.annotate(location.address, xy=pos, xycoords='data',\n",
    "             xytext=(0.25, 0.5), textcoords='figure fraction',color=\"red\",\n",
    "             arrowprops=dict(arrowstyle=\"->\"))\n",
    "ax.set_aspect('equal')\n",
    "# Set a point\n",
    "circ = Circle((pos[0],pos[1]),60,color=\"red\")\n",
    "ax.add_patch(circ)\n",
    "fig.set_size_inches(18.5, 10.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now, well, OSM is hard.\n",
    "* This is I wanted to show where to programaticaly query OSM for all sorts of data.\n",
    "* Turns out that it is a lot harder than it should be, especially if you want to work with GeoJson.\n",
    "* Out of scope for this talk, but let's punt.\n",
    "* OSM has a feature called Overpass. It is like the most convoluted Google maps ever using a very complex query language that I still don't grok. \n",
    "* We're going to use it to get all of the parks in Portland as GeoJSON using the web interface called Overpass Turbo.\n",
    "* [Let's take a look at that.](https://overpass-turbo.eu/)\n",
    "* Here's the query to run. Then export as GeoJSON\n",
    "```\n",
    "[bbox:{{bbox}}][timeout:1800];\n",
    "way[\"leisure\"=\"park\"];map_to_area->.a;\n",
    "foreach(\n",
    "  (._;>;);\n",
    "  is_in;\n",
    "  way(pivot)[\"leisure\"=\"park\"];\n",
    "  out geom;\n",
    ");\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's load up our park data\n",
    "* Load the file using GeoPandas (some syntactic sugar on Pandas). \n",
    "* Also load the raw json, and chunk out each park.\n",
    "* Update the area value because there is no value, not really useful except as a proxy measurement.\n",
    "* Update and sort our data frame. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "fname = \"./portland_parks_small.geojson\"\n",
    "park_df = gpd.read_file(fname)\n",
    "portland_parks = json.load(open(fname,'r'))\n",
    "# raw geojson works better with GDAL\n",
    "geojson = [p for p in portland_parks[\"features\"]]\n",
    "# no area out of the box\n",
    "p = [p.area for p in park_df[\"geometry\"].tolist()]\n",
    "park_df[\"area\"] = pd.Series(p)\n",
    "park_df[\"geojson\"] = pd.Series(geojson)\n",
    "park_df.sort_values(['area',], ascending=[1])\n",
    "park_df.head()\n",
    "#len(park_df)\n",
    "#print park_df[\"wikipedia\"].dropna()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now we'll update our slippy map.\n",
    "* Just toss the aois in, just like our scene footprints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for p in portland_parks[\"features\"]:\n",
    "    feat = {'geometry':p[\"geometry\"],\"properties\":{\n",
    "            'style':{'color': \"#00FF00\",'fillColor': \"#00FF00\",'fillOpacity': 0.0,'weight': 1}},\n",
    "            'type':u\"Feature\"}\n",
    "    # convert to geojson\n",
    "    gjson = GeoJSON(data=feat)\n",
    "    # add it our map\n",
    "    m.add_layer(gjson)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now let's find the big parks.\n",
    "* The pandas dataframe can have multiple enteries per park.\n",
    "* We can use the group by command to sum up these disparate areas. \n",
    "* Finally we'll output the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "park_sz = park_df.groupby(\"name\").sum() \n",
    "park_sz = park_sz.sort_values(by='area',ascending=[0])\n",
    "display(park_sz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now to the meat of the problem.\n",
    "* Our goal is to get each park as a small image so we can analyze it.\n",
    "* We'll write a function to create a geojson file from our big geojson file\n",
    "* We'll also write a function that takes in our scene list, an input geojson file, and calls gdal warp to generate our small park image. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def write_geojson_by_name(df,name,outfile):\n",
    "    \"\"\"\n",
    "    Take in a dataframe, a park name, and an output file name \n",
    "    Save the park's geojson to the specified file.\n",
    "    \"\"\"\n",
    "    temp = df[df[\"name\"]==name]\n",
    "    to_write = {\"type\": \"FeatureCollection\",\n",
    "                \"features\": temp[\"geojson\"].tolist()}\n",
    "    with open(outfile,'w') as fp:\n",
    "        fp.write(json.dumps(to_write))\n",
    "def crop_scenes_to_geojson(geojson,scenes,out_name):\n",
    "    \"\"\"\n",
    "    Take in a geojson file, a list of scenes, and an output name\n",
    "    Call gdal and warp the scenes to match the geojson file and save the results to outname.\n",
    "    \"\"\"\n",
    "    commands = [\"gdalwarp\", # t\n",
    "           \"-t_srs\",\"EPSG:3857\",\n",
    "           \"-cutline\",geojson,\n",
    "           \"-crop_to_cutline\",\n",
    "           \"-tap\",\n",
    "            \"-tr\", \"3\", \"3\"\n",
    "           \"-overwrite\"]\n",
    "    for tiff in scenes:\n",
    "        commands.append(tiff)\n",
    "    commands.append(out_name)\n",
    "    print \" \".join(commands)\n",
    "    subprocess.call(commands)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's put it all together\n",
    "* We're going to use the scenes we downloaded earlier as our input and build a little image for every park in Portland!\n",
    "* We just have to make a few file names and call the functions above.\n",
    "* If we really wanted to get fancy we could do this for every image that has our park and make a sick movie or lots of different types of a analysis. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_json_files = []\n",
    "tif_file_names = []\n",
    "unique_park_names = list(set(park_df[\"name\"].tolist()))\n",
    "for name in list(unique_park_names):\n",
    "    # Generate our file names \n",
    "    geojson_name = \"./parks/\"+name.replace(\" \",\"_\")+\".geojson\"\n",
    "    tif_name = \"./parks/\"+name.replace(\" \",\"_\")+\".tif\"\n",
    "    # write geojson\n",
    "    write_geojson_by_name(park_df,name,geojson_name)\n",
    "    # write to park file\n",
    "    crop_scenes_to_geojson(geojson_name,tiff_files,tif_name)\n",
    "    # Save the results to lists\n",
    "    geo_json_files.append(geojson_name)\n",
    "    tif_file_names.append(tif_name)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's take a look at the first few parks!\n",
    "* matplotlib and tifs can be a bit heavy handed. \n",
    "* imma teach you a protip use [image magick](https://www.imagemagick.org/script/index.php) and the built in image display.\n",
    "* Use subprocess to tell imagemagick to convert tifs to jpg.\n",
    "* Then load and display the images. \n",
    "* WARNING: do not use imagemagick to modify geotiffs!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "magic = [\"mogrify\",\"-format\", \"jpg\", \"./parks/*.tif\"]\n",
    "subprocess.call(magic)\n",
    "for p in tif_file_names[0:20]:\n",
    "    print p\n",
    "    display(Image(p.replace('tif','jpg')))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Not let's do some quick analytics -- your code goes here.\n",
    "* For completeness let's do some basic image processing.\n",
    "* For each of parks we are going to calculate the average \"greeness\" per pixel over the other two channels.\n",
    "* We do this as it controls for white pixels, like clouds.\n",
    "* Since there are black pixels we'll have to controll for that by only using only the non-black pixels. \n",
    "* We'll use numpy here, but scikit image and OpenCV have many more features. \n",
    "* It is also worth noting that the visual product is probably only useful for calculating areas. If you want to do real science use the Analytics products.\n",
    "* The real way to do this is to calculate a [Normalized Difference Vegetation Index (NDVI)](https://en.wikipedia.org/wiki/Normalized_Difference_Vegetation_Index) using the analytic product.\n",
    "* [Here is an example of NDVI calculations](https://github.com/planetlabs/notebooks/blob/master/ndvi/ndvi_planetscope.ipynb). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def load_image3(filename):\n",
    "    \"\"\"Return a 3D (r, g, b) numpy array with the data in the specified TIFF filename.\"\"\"\n",
    "    path = os.path.abspath(os.path.join('./', filename))\n",
    "    if os.path.exists(path):\n",
    "        with rasterio.open(path) as src:\n",
    "            b,g,r,mask = src.read()\n",
    "            return np.dstack([b, g, r])\n",
    "        \n",
    "\n",
    "def get_avg_greeness(filename):\n",
    "    retVal = -1.0\n",
    "    try:\n",
    "        # load the image\n",
    "        img = load_image3(filename)\n",
    "        if img is not None:\n",
    "            # add all the channels together, black pixels will still be zero\n",
    "            # this isn't a perfect method but there are very few truly black spots \n",
    "            # on eart\n",
    "            black_like_my_soul = np.add(np.add(img[:,:,0],img[:,:,1]),img[:,:,2])\n",
    "            # sum up the not black pixels\n",
    "            not_black = np.count_nonzero(black_like_my_soul)\n",
    "            # sum up all the green\n",
    "            img = np.array(img,dtype='int16')\n",
    "            total_green = np.sum(img[:,:,1]-((np.add(img[:,:,0],img[:,:,2])/2)))\n",
    "            # calculate our metric\n",
    "            if total_green != 0 and not_black > 0:\n",
    "                retVal = total_green / float(not_black)\n",
    "        return retVal\n",
    "    except Exception as e:\n",
    "        print e\n",
    "        return -1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "greens = [get_avg_greeness(f) for f in tif_file_names]\n",
    "print greens"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "paired = zip(tif_file_names,greens)\n",
    "paired.sort(key=(lambda tup: tup[1]))\n",
    "paired.reverse()\n",
    "labels = [p[0][8:-4].replace(\"_\",\" \") for p in paired]\n",
    "data = [p[1] for p in paired]\n",
    "plt.figure(figsize=(20,6))\n",
    "xlocations = np.array(range(len(paired)))+0.5\n",
    "width = 1\n",
    "plt.bar(xlocations, data, width=width)\n",
    "plt.yticks(range(-1,25,1))\n",
    "plt.xticks(xlocations+ width/2, labels)\n",
    "plt.xlim(0, xlocations[-1]+width*2)\n",
    "plt.ylim(-2,np.max(data)+1)\n",
    "plt.title(\"Greeness over Average Red and Blue Per Park\")\n",
    "plt.gca().get_xaxis().tick_bottom()\n",
    "plt.gca().get_yaxis().tick_left()\n",
    "xa = plt.gca()\n",
    "xa.set_xticklabels(xa.xaxis.get_majorticklabels(), rotation=90)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's take a look at what this looks like in terms of images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imgs = [p[0] for p in paired]\n",
    "for p in imgs[0:35]:\n",
    "    print p[8:-4].replace(\"_\",\" \")\n",
    "    display(Image(p.replace('tif','jpg')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now let's plot over our slippy map. \n",
    "* We'll calculate a non-linear opacity per park and then use that for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opacity_map = {}\n",
    "gmax = np.max(greens)\n",
    "gmin = np.min(greens)\n",
    "# this is a nonlinear mapping\n",
    "opacity = [np.clip((float(g**2)-gmin)/float(gmax-gmin),0,1) for g in greens]\n",
    "for op,name in zip(opacity,imgs):\n",
    "    opacity_map[name]=op\n",
    "\n",
    "m = Map(\n",
    "    center=center, \n",
    "    zoom=zoom,\n",
    "    default_tiles = planetMapTiles # Uncomment to use Planet.com basemap\n",
    ")\n",
    "dc = DrawControl(\n",
    "    polygon = polygon,\n",
    "    rectangle = rectangle\n",
    ")\n",
    "# Initialize an action counter variable\n",
    "actionCount = 0\n",
    "AOIs = {}\n",
    "\n",
    "# Register the draw controls handler\n",
    "def handle_draw(self, action, geo_json):\n",
    "    # Increment the action counter\n",
    "    global actionCount\n",
    "    actionCount += 1\n",
    "    # Remove the `style` property from the GeoJSON\n",
    "    geo_json['properties'] = {}\n",
    "    # Convert geo_json output to a string and prettify (indent & replace ' with \")\n",
    "    geojsonStr = json.dumps(geo_json, indent=2).replace(\"'\", '\"')\n",
    "    AOIs[actionCount] = json.loads(geojsonStr)\n",
    "    \n",
    "# Attach the draw handler to the draw controls `on_draw` event\n",
    "dc.on_draw(handle_draw)\n",
    "m.add_control(dc)\n",
    "m\n",
    "\n",
    "for p in portland_parks[\"features\"]:\n",
    "    t = \"./parks/\"+p[\"properties\"][\"name\"].replace(\" \",\"_\") + \".tif\"\n",
    "    feat = {'geometry':p[\"geometry\"],\"properties\":{\n",
    "            'style':{'color': \"#00FF00\",'fillColor': \"#00FF00\",'fillOpacity': opacity_map[t],'weight': 1}},\n",
    "            'type':u\"Feature\"}\n",
    "    # convert to geojson\n",
    "    gjson = GeoJSON(data=feat)\n",
    "    # add it our map\n",
    "    m.add_layer(gjson)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
